1 Setup

# load libraries
library(tidyverse) # dplyr, tidyr, ggplot
library(openxlsx) # reading and writing excel file
library(knitr) # generating reports
library(plotly) # interactive plots
library(CHNOSZ) # stability diagrams
data(thermo)
opts_chunk$set(
  collapse = TRUE, comment = "#>",
  dev=c("png", "pdf"), dev.args=list(pdf = list(encoding="WinAnsi", useDingbats=FALSE))
  )

2 Mineral stability diagram K\(_2\)O - Al\(_2\)O\(_3\) - SiO\(_2\) - H\(_2\)O system

par(mfrow = c(2, 2))
res <- 200
fill <- "terrain"

## K2O-Al2O3-SiO2-H2O, 25 degree C, 1 bar
## Steinmann et al., 1994 (http://ccm.geoscienceworld.org/content/42/2/197)
## Garrels and Christ, p. 361 (http://www.worldcat.org/oclc/517586)
## https://wwwbrr.cr.usgs.gov/projects/GWC_coupled/phreeqc/html/final-75.html
basis(c("Al+3", "pseudo-H4SiO4", "K+", "H2O", "H+", "O2"))
species(c("gibbsite", "muscovite", "kaolinite", "pyrophyllite", "K-feldspar"))
a <- affinity(H4SiO4 = c(-6, -2, res), `K+` = c(-3, 6, res))
#> energy.args: temperature is 25 C
#> energy.args: pressure is Psat
#> energy.args: variable 1 is log_a(H4SiO4) at 200 values from -6 to -2
#> energy.args: variable 2 is log_a(K+) at 200 values from -3 to 6
#> subcrt: 11 species at 298.15 K and 1 bar (wet)
diagram(a, ylab = ratlab("K+"), fill = fill, yline = 1.7)
#> balance: from moles of Al+3 in formation reactions
#> diagram: plotting A/(2.303RT) / n.balance (maximum affinity method for 2-D diagrams)
title(main = syslab(c("K2O", "Al2O3", "SiO2", "H2O")))
#> Warning in title(main = syslab(c("K2O", "Al2O3", "SiO2", "H2O"))): font
#> metrics unknown for Unicode character U+2013

#> Warning in title(main = syslab(c("K2O", "Al2O3", "SiO2", "H2O"))): font
#> metrics unknown for Unicode character U+2013

#> Warning in title(main = syslab(c("K2O", "Al2O3", "SiO2", "H2O"))): font
#> metrics unknown for Unicode character U+2013

#> Warning in title(main = syslab(c("K2O", "Al2O3", "SiO2", "H2O"))): font
#> metrics unknown for Unicode character U+2013
#> Warning in title(main = syslab(c("K2O", "Al2O3", "SiO2", "H2O"))):
#> conversion failure on '–' in 'mbcsToSbcs': dot substituted for <e2>
#> Warning in title(main = syslab(c("K2O", "Al2O3", "SiO2", "H2O"))):
#> conversion failure on '–' in 'mbcsToSbcs': dot substituted for <80>
#> Warning in title(main = syslab(c("K2O", "Al2O3", "SiO2", "H2O"))):
#> conversion failure on '–' in 'mbcsToSbcs': dot substituted for <93>
#> Warning in title(main = syslab(c("K2O", "Al2O3", "SiO2", "H2O"))): font
#> metrics unknown for Unicode character U+2013
#> Warning in title(main = syslab(c("K2O", "Al2O3", "SiO2", "H2O"))):
#> conversion failure on '–' in 'mbcsToSbcs': dot substituted for <e2>
#> Warning in title(main = syslab(c("K2O", "Al2O3", "SiO2", "H2O"))):
#> conversion failure on '–' in 'mbcsToSbcs': dot substituted for <80>
#> Warning in title(main = syslab(c("K2O", "Al2O3", "SiO2", "H2O"))):
#> conversion failure on '–' in 'mbcsToSbcs': dot substituted for <93>
#> Warning in title(main = syslab(c("K2O", "Al2O3", "SiO2", "H2O"))): font
#> metrics unknown for Unicode character U+2013
#> Warning in title(main = syslab(c("K2O", "Al2O3", "SiO2", "H2O"))):
#> conversion failure on '–' in 'mbcsToSbcs': dot substituted for <e2>
#> Warning in title(main = syslab(c("K2O", "Al2O3", "SiO2", "H2O"))):
#> conversion failure on '–' in 'mbcsToSbcs': dot substituted for <80>
#> Warning in title(main = syslab(c("K2O", "Al2O3", "SiO2", "H2O"))):
#> conversion failure on '–' in 'mbcsToSbcs': dot substituted for <93>
legend("bottomleft", describe.property(c("T", "P"), c(25, 1)), bty = "n")

3 K feldspar weathering to gibbsite

KAlSi\(_3\)O\(_8\) + H\(^+\) + 7H\(_2\)O \(\leftrightharpoons\) Al(OH)\(_3\) + 2K\(^+\) + 3H\(_4\)SiO\(_4\)

3.1 Initial conditions

mol_gibbs_t0 <- 0 # moles gibbsite at t0
mol_l_H_t0 <- 1e-4 # moles per liter H+ t0
mol_l_K_t0 <- 1e-6 # moles per liter K+ t0
mol_l_H4SiO4_t0 <- 1e-6 # moles per liter H4SiO4 t0 
#mol_kaol_t0 <- 0 # moles kaolinite t0

rxn_prog_var <- 7.47e-7 # reaction progress variable is proportional to the amount of reactant mineral dissolved and defines the extent of reaction. Selection of reaction progress step is arbitrary. units are moles per kg of solution

#coefficients of change for weathering of k spar to gibbsite (equation above)
co_spar_1 <- -1
co_H_1 <- -1
co_H2O_1 <- -7
co_gibbs_1 <- 1
co_K_1 <- 1
co_H4SiO4_1 <- 3
steps <- seq(0, 26, by = 0.01)

k_spar_to_gibbs <- as.data.frame(steps)

#k_spar_weathering <- k_spar_weathering %>% rename("f_rem" = "f_rem_list")
k_spar_to_gibbs <- k_spar_to_gibbs %>% mutate(mol_l_H = mol_l_H_t0 + co_H_1 * steps * rxn_prog_var, mol_l_K = mol_l_K_t0 + co_K_1 * steps * rxn_prog_var, mol_l_H4SiO4 = mol_l_H4SiO4_t0 + co_H4SiO4_1 * steps * rxn_prog_var, mol_gibbs = mol_gibbs_t0 + co_gibbs_1 * steps * rxn_prog_var, log_H4SiO4 = log10(mol_l_H4SiO4), log_K_H_ratio = log10(mol_l_K/mol_l_H))
plot_k_spar_to_gibbs <- k_spar_to_gibbs %>% ggplot(aes(x=log_H4SiO4, y = log_K_H_ratio, color = steps))+
  geom_point() +
  scale_x_continuous(limits = c(-6, -2))+
  scale_y_continuous(limits = c(-3, 6))

plot_k_spar_to_gibbs %>% ggplotly()

4 K feldspar weathering

2KAlSi\(_3\)O\(_8\) + 2H\(^+\) + 9H\(_2\)O \(\leftrightharpoons\) Al\(_2\)Si\(_2\)O\(_5\)(OH)\(_4\) + 2K\(^+\) + 4H\(_4\)SiO\(_4\)